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Abstract 

A novel methodology for accelerating the solution of PDE-constrained optimization is introduced. 
It is based on an offline construction of database of local ROMs and an online interpolation within 
the database. The online flexibility of the ROM database approach makes it amenable to speeding- 
up optimization-intensive applications such as robust optimization, multi-objectives optimization, 
and multi-start strategies for locating global optima. The accuracy of the ROM database model can 
be tuned in the offline phase where the database of local ROMs is constructed through a greedy 
procedure. In this work, a novel greedy algorithm based on saturation assumption is introduced to 
speed-up the ROM database construction procedure. The ROM database approach is applied to a 
realistic wing design problems and leads to a large online speed-up. 
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1 Introduction 

Many physical and social phenomena can be described by Partial Differential Equations (PDEs) and ac¬ 
curately simulated thanks to advances in numerical analysis and computer technology. PDE-constrained 
optimization problems arise in numerous important applications: design optimization, inverse problems, 
and optimal control. In spite of the importance of PDE-constrained optimization, many difficulties are 
known for the process of solving it. First of all, it is hard to find a robust PDE solver that works for 
dramatic changes in parameter values. This issue, although it has been resolved to a certain extent due 
to active research on this topic, remains an ongoing research topic as the complexity of applications 
requiring PDE modeling increases. Second, a PDE solve can be very expensive for complex problems. 
This causes the optimization process to be impractically long due to the many queries to the PDE solver 
involved. This second difficulty can be resolved by replacing the PDE with a surrogate model such as 
a projection-based Reduced Order Model (ROM) that lies within the subspace spanned by a Reduced 
Order Basis (ROB) [34, 29]. Unfortunately, it has been shown that ROMs that are constructed for a 
given value of parameters are in general not robust with respect to parameter changes [15]. 

In the context of optimization, the robustness of a ROM can be addressed in three different ways. 
First, a global ROM that is globally accurate in the parameter space can be constructed offline and 
used in the optimization process online [8, 37, 26]. Because there is no call to the PDE solver in the 
online phase, the optimization process can be accelerated tremendously. However, the accuracy of the 
global ROM over the whole parameter space depends heavily on the size of the ROB. The bigger the 
parameter space is, the larger the size of ROB must be to have a sufficient accuracy on the parameter 
space. The second way of improving the robustness of a ROM in optimization is to adaptively update 
the reduced-order basis of the global ROM [38, 39, 40, 16]. As the optimization progresses, some basis 
vectors in the ROB are eliminated and some new vectors are added in order to improve the accuracy 
of the global ROM around the current point. Due to the locality of the global ROM, this progressive 
approach gives a better accuracy than using one global ROM over the whole parameter space. However, 
the optimization process becomes slow because updating the global ROM adaptively requires calling the 
PDE solves. 

The third way of improving the robustness of a ROM in optimization is to use a database of local 
parameterized ROMs and interpolate the ROMs to quickly generate the new ROMs at non-populated 
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data points [3, 1]. The attractiveness of this approach is two-fold: the availability of local ROMs with a 
small size of ROB over the whole parameter space and the avoidance for constructing a new ROM within 
the optimization process. Additionally, the ROM database approach gives great online flexibility. For 
example, the ROM database approach is appealing when the multiple optimization solves are necessary 
such as robust optimization, multi-objectives optimization, and global optimization with multiple starts. 
It is because the database of local ROMs can be reused without additional offline phase. In the context 
of multidisciplinary problems, different databases from different disciplinary can also be easily combined. 
Finally, the fast online phase in the ROM database approach enables real-time optimization for time- 
critical applications. 

This paper develops a methodology for a ROM database approach in the optimization context. 
A novel, cost-efficient greedy algorithm for constructing a database is introduced and compared with 
existing greedy algorithms. The new greedy algorithm can be applicable not only to the construction 
of local ROM database, but also to constructing a global ROM where snapshots from parameter space 
need to be chosen in a smart way. It also explains how to obtain ROMs and their sensitivities at non 
populated data points. Finally, the proposed database approach is applied to a realistic wing design 
problems. 

The rest of the paper is organized as follows. In Section 2, the optimization problem of interest 
and projection-based model reduction are presented. The main idea of the proposed methodology for 
solving a PDF-constrained optimization is described in Section 3. A new approach for constructing a 
database VB is presented and compared with existing approaches in Section 4. Section 5 presents the 
consistent interpolation of ROMs on matrix manifolds and derives the gradients of the ROM interpolant 
with respect to design parameters fi. Section 6 demonstrates the solution procedure for the ROM 
database model-constrained optimization in an aeroelastic wing shape optimization. Section 7 concludes 
the paper. 

2 Problem statement 

The following optimization problem is considered: 


minimize 



subject to 

c(w(/i),/i) < 0 


( 1 ) 


where /(•,•) G M is an objective function, €(•,•) G defines Nc inequality constraints, fji e V G 
is a vector of optimization variables, P is a parameter space, and w G is a vector of state 
variables solution of a parameterized linear system: 

( 2 ) 

This linear system typically arises from the linearization of a non-linear PDF around a nominal condi¬ 
tion. The optimization problem (1) only handles the parameter vector (jl and not the vector of state 
variables w. This problem therefore pertains to the framework of Nested Analysis and Design (NAND). 
Alternative framework for solving a PDF-constrained optimization is Simultaneous Analysis and Design 
(SAND) [8, 11, 13]. In a gradient-based optimization algorithm, the first derivatives of the objective 
function /(w(/L6), /l^) and each of the constraints fi) for i = 1, • • • , A'c need to be computed. In 

general, if q denotes a generic quantity of interest whose derivative with respect to is required such as 
/ and Q, the chain rule leads to 

^(wW.M) = (3) 

Fquation (2) are also differentiated for each parameter /r^, i = 1 • • • , leading to a linear system: 
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Substituting the sensitivity solution of (4) into (3) leads to 


(w(/Lfc),/Lfc) 


w(ax),Ax) + A{fi) 


1 f dh 




There are two approaches for computing the set of sensitivities {^(w(/L6), 

1. In the direct approach^ the state sensitivities are first computed by solving the linear system 

(5), then the sensitivities ^ can be evaluated by (3). 

2. In the adjoint approach^ the adjoint vector Xq{/J^) is first computed by solving the following linear 
system A{/jl)^X q{/ji) = ^(w(/L6),/l^) and all the sensitivities ^ are compute as 

(f-*"’ - Is*"***"’) . i = 1. ■ ■ ■ (6) 


The direct approach requires linear solves while the adjoint approach requires A^c + 1 solutions. 
Hence, if < 1+A^c the direct approach is preferable while the adjoint approach is preferable otherwise. 
Problem (1) can then be solved by a gradient-based nonlinear optimization algorithm such as Sequential 
Quadratic Programming (SQP) [19] together with a quasi-Newton approximation of the Hessian matrix, 
the trust-region method [14], or the interior-point method [36]. 

The solution of the linearized PDE and its associated sensitivity or adjoint equations is compu¬ 
tationally expensive as it requires the solution of linear systems of dimension N^. To alleviate that 
cost, projection-based model reduction reduces the dimension of the system to be solved by reducing 
the dimensionality of the state w(/L 6). For that purpose, a pre-computed reduced-order basis (ROB) 
y{pi) G spanning a /c-dimensional subspace S{pl) C is defined and the state approximated 

as: 

w{^l) (7) 

where k <C and Wr(/Lt) € denotes the reduced coordinates of the state w(/Lt) in terms of the ROB 
y{pi). The state approximation (7) introduces a (usually) non-zero residual r(w^,/L6) associated with 
the parameterized linear system (2) defined as: 

r(w^,/L6) = A{pL)\{pL)Wr - h{pL). (8) 

This residual is then enforced to be orthogonal to a second ROB, W(/L6) G as W(/L6)^r(w^, fi) = 0 

resulting in the reduced linear system of dimension k: 

A^(/Lfc)w^(/Lfc) = W(/Lfc), (9) 

where the parameterized reduced order matrix and the reduced vector are defined as A^(/L6) = 
WA{pl)\( fi) and hr {pi) = W(/L6)^b(/L6), respectively. This results in a Petrov-Galerkin projection. 
If W(/L6) = V(/L6), this is a Galerkin projection. The optimization problem (1) can then be replaced by 
the following cheaper problem involving the reduced coordinates only: 

( 10 ) 

where Wr{pi) is a solution of the reduced linear system of equtions, Ar{pi)wr{pi) = hr{pi). Gonstructing 
a set of reduced bases (V(/L6), W(/L6)) for a given parameter pi can be done by Proper Orthogonal 
Decomposition (POD) [34], Balanced Truncation [29] or Moment Matching [18]. However, all of these 
approaches involve intensive computations in order to construct (V{pi),W{pi)). To address this issue, 
ROB and ROM interpolation approaches have been developed in [3, 1]. All of these approaches proceed 
by pre-computing a database of reduced operators 

PB={(Me,V(Rc),W(Rc))}t^ or VB = {{fX,, . (11) 
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Offline 


Online 



Figure 1: Methodology flow chart for the offline and online phases 


The set VB is then interpolated to cheaply construct reduced order bases and W{iJ,'^) or reduced 

order models and for a parameter ^V. A database of RGBs {{V{/^i){/^i))}c=i 

can be interpolated on the Grassmannian manifold [2]. However, the interpolation of the RGBs is more 
expensive than the interpolation of the RGMs because the interpolation on the Grassmannian manifold 
requires Singular Value Decompositions (SVDs) of matrices scaling with the size N^. Therefore, this 
paper focuses on the inexpensive interpolation of the RGMs involving reduced operators only. 

In order to solve PDE-constrained Gptimization (1) efficiently with sufficient accuracy, the RGM- 
constrained optimization problem (10) can be solved, instead, using the database approach where a RGM 
interpolation strategy is used for robustness and efficiency. Section 3 presents the proposed methodology 
of solving the RGM-constrained optimization (10). 


3 ROM-constrained optimization 

The proposed procedure for solving the RGM-constrained optimization problem (10) can be divided 
into two separate phases: an offline phase followed by an online phase. In the offline phase, the RGM- 
database VB is constructed. In the online phase, the RGM-constrained optimization problem (10) is 
then solved by a gradient-based optimization algorithm where the function values and its derivatives are 
computed by RGM interpolation of the elements in the database VB. Figure 1 schematically describes 
the methodology flow chart both for the offline and online phases. In the figure, the offline phase chooses 
six points in VB. At each point in VB^ the corresponding local RGMs are stored (see Eq. (11)). These 
RGMs are then interpolated in the online phase in order to obtain RGMs at fi* ^ VB as in Eigure 1. 
The interpolant RGMs at is then used to compute and pass and |J(m*) to a gradient-based 

optimizer and the optimizer iterates until it converges. 

The efficiency and accuracy of the methodology depends on how well the database is constructed in 
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the offline phase. For example, if a large number of points of V are included in PS, the accuracy of the 
model raises up to the level of the local ROM accuracy, but the efflciency is small as the offline phase 
is extremely expensive. On the other hand, if there are very few points in PS, then both the offline 
and online phases are fast, but the accuracy of the ROM is low. These two extreme cases illustrate 
the importance of an “optimal” database construction. The “optimal” database PS* can be abstractly 
defined as the one that maximizes the product of efflciency and accuracy: 

PS* = argmax efflciency (PS) • accuracy (PS). (12) 

VB 

Finding the optimal database PS* is a challenging task in the offline phase. 

4 Database construction 

There are two main approaches for constructing a database: by a priori sampling approach and by 
adaptive sampling. The a priori sampling approach requires no information about error or accuracy of 
the model on P whereas the adaptive sampling approach requires knowledge about error of the model. 
The a priori sampling approach tries to select samples in PS so that PS is a good representation of 
P topologically or statistically. Examples of a priori sampling include full factorial design and latin 
hypercube sampling. On the other hand, the adaptive sampling approach updates PS in a way that the 
update produces the “maximum” increase in accuracy of the model. Greedy algorithms are widely used 
in the adaptive sampling approach [8, 9, 10, 21, 24, 31, 35]. Because the a priori sampling approach does 
not depend on the model, it leads to a fast construction of the database but tends to include unnecessary 
samples. On the other hand, the adaptive sampling approach tends to construct a database that is closer 
to an optimal database. However, the procedure of adaptive sampling is more computationally expensive 
than the one of a priori sampling approach because the adaptive sampling approach relies on repeated 
evaluations of the errors of the model (or some error indicators) to assess the accuracy of the model. 

Section 4.1 briefly summarizes the full factorial design and latin hypercube sampling approaches. 
Section 4.2 describes the proposed adaptive sampling technique based on a greedy algorithm. Both 
sections focus on aspects of the sampling methods in the context of the ROM database model. 

4.1 A priori sampling 

In a full factorial design, each variable domain is divided into several factor levels. Every combination 
of the factor levels is then included in the database VB. An example of full factorial design in is 
depicted in Eigure 2 (left) where all the factor levels are set to five. By the nature of the design, all 
the points in VB produced by full factorial design are uniformly spaced, which is a desired property 
in the context of interpolation because interpolation of points with irregular spacing may result in an 
ill-conditioned system of equations. However, the number of points in VB increases exponentially as 
the size of the parameter space increases. Eor example, if P C and each variable’s factor level is 
three, then 729 points are in VB which becomes expensive. Additionally, it is hard to determine a priori 
the appropriate size of the database for a certain accuracy of the model because the full factorial design 
does not use any information about the accuracy of the model. 

Latin hypercube sampling tries to overcome the issue of oversampling in VB by randomly generating 
points in a way that the generated points are well distributed in VB. In latin hypercube sampling, each 
variable in P C is uniformly divided into a same number of partitions (e.g., M partitions for each 
variable). It then enforces to have one sample point in each dimension (see Eigure 2 (right)). Note that 
the number of points in VB only depends on the number of partitions, M, not on N^. This is a desirable 
property in order to avoid the curse of dimensionality. However, the accuracy of the model now strongly 
depends on M. The higher the M value is, the more accurate the model is. Thus, latin hypercube 
sampling does not provide a complete freedom from the curse of dimensionality because the higher 
is, the higher M is required for the accuracy of the model. Additionally, the points in VB generated by 
latin hypercube sampling are not guaranteed to have even spacing. This may cause issues in accuracy 
of the ROM database approach away from sample points. Also, latin hypercube sampling is unlikely 
to choose corner points, so some points in V have to be evaluated by extrapolation only. Einally, it is 
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Figure 2: Two a priori sampling approaches: full factorial design (left), latin hypercube sampling (right). 


hard to know a priori the appropriate M value for a certain accuracy of the model as in the case of full 
factorial design as latin hypercube sampling does not use any information about the accuracy of the 
model. 

4.2 Adaptive sampling 

4.2.1 Classical greedy procedure 

Adaptive sampling approaches make use of error estimates associated with the ROM to populate the 
database VB iteratively. Greedy algorithms are widely used as simple adaptive sampling approach. At 
each iteration, the greedy algorithm selects a point where the maximum error occurs in V and reduces the 
error by including the point in VB. More formally, denoting by VBn a ROM database with Np samples, 
let q{ii]VBn ) denote an error indicator for the VBn -based ROM-interpolation model dit fj, eV. Let 
S also define a candidate set of Ns test points in V. The candidate set S must include enough points 
to represent V well. One can either generate E by full factorial design with a large value of factor level 
for each variable or by latin hypercube sampling with a large value of M. The greedy algorithm first 
picks a point fi^ in E and builds the initial database VBi = At iteration Ap, the 

algorithm computes ) for i = 1,..., Ns and adds the maximum error point in E to VBn to 

form VBnp+ 1 - The greedy algorithm repeats this process until the maximum error indicator among the 
points in E is less than a threshold Ctoi • 

In practice, there are two types of error indicators: 1) error bounds A(/L6) and 2) indicators based on 
the norm of the residual ||r(w)||. Error bounds A(w) rigorously satisfy 

||w* — Vw^ll < A(Vw^), Vw^ G (13) 

where w* denotes the PDE solution. Such error bounds have been derived for elliptic [32, 35], parabolic 
[21] and hyperbolic [22] PDEs, linear time invariant systems [23], and nonlinear eigenvalue problems 
[12]. However, these error bounds typically require the computation of an inf-sup constant and are not 
usually tight [5]. They are therefore currently limited to the aforementioned classes of equations. The 
norm of the residual is another popular alternative to error bounds [9, 10] in cases where such an 
error bound does not exist or is not a tight indicator of the error. 

The requirement of computing q{ii]VBnp) for aW fj, eE at each iteration Np can be very expensive 
when the number of candidates Ns is large. Two recent papers address this issue [31, 24]. The authors 
in [31] propose to use a surrogate model of error surface to find the point of the global maximum error 
indicator in V. The authors in [24] propose an improved greedy algorithm by adopting the concept of the 
saturation assumption, which is widely used in the adaptive finite element mesh refinement algorithm. 
Alternative adaptive greedy algorithms are also proposed in the subsequent section. 


6 














































Algorithm 1 Saturation assumption-based adaptive greedy algorithm for the ROM-interpolation model 
Input: A candidate set S C P, the candidate set size As, marginal factor 7, a tolerance etoi > 0, an 
initial saturation constant 77 , the subset size An 
Output: A ROM database VB^p 

1: Choose an initial parameter value (jl^ G S, compute Ap{/j.-^^) and hp{/j,-^) 

2: Set VBi = {{fXj, 

3: Set ^profiiKM) = £'proflie(M) = oo for all G S and = oo 

4: while > €tol do 

5: Generate a random subset 11 of S 

6: Set = 1 and rtemp(j) = 0, j = 1,..., A^n 

7: for G njVp,i = do 

8: if TsQpromeiBj) > Qmayi and TgQpromeiBj) > ^tol then 

9: Set = gpromeiBj) 

10: Compute g(pj-, VBn ^), set Pproflie (Bj) = eiPj', ) 

11- if £’proflk(Mi) 7^ 00 and PprofiKMi) > etoi then 

12: Set TtempO') = max(l, 

13: end if 

14: if p(/Ltj;X>BjVp) > Pmax then 

15: Set = gifipVBNj and Pn^+i = fij 

16: end if 

17: end if 

18: end for 

19: if fi’max < etoZ then 

20: Perform a sanity check 

21: end if 

22: if max(Ttemp) ^ 1 then 

23: Set Ts = max(Ttemp) 

24: end if 

25: Compute and 

26: Set VBnp + 1 = 'B>Bnp U {(^ATp + l? Ar(p-Arp + l): W(MArp + l))} 

27: Np i — Np 1 

28: end while 


4.2.2 Random greedy procedure 

In order to speed up the process of the greedy algorithm, a random subset of E can be evaluated at 
each iteration, instead of the whole candidate set E. Let 11 at denote a subset of E with size An ^ A^ 
randomly selected from E at Iteration Np. The subset Hat is updated in each greedy iteration in order 
to explore the parameter space V. Additionally, only the points away from VBnp are included in IlATp. 
This condition can be reasoned from the fact that the true maximum error is likely to happen at the 
point away from VBnp- The process is repeated until a convergence condition (i.e., the maximum error 
estimate p^ax than a desirable convergence tolerance e^oz) is satisfied. If the convergence condition 

is satisfied at the end of the current greedy iteration, a sanity check is done by evaluating error indicators 
at points in a random subset IIaTp of larger size. The sanity check is required to ensure that the database 
VBnp does not have any important missing points in V. If the convergence condition < ctoi is still 
satisfied after the sanity check, then it returns VBnp- If not, the point that gives the maximum error 
indicator in the sanity check is added to the database and the greedy procedure is continued. 

4.2.3 Saturation constant assumption-based greedy procedure 

The random greedy procedure can be further accelerated by applying the saturation assumption- 
based filtering proposed by [24]. It is detailed as follows. 
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Definition 1. Saturation Constant 

Let Q^fi] VBn ) denote an error indieator depending on a parameter /j, and a database VBn with nested 
databases satisfying VBn C VBm for all 1 < Np < Mp. The Saturation Constant Ts > 0 is defined 
as 


g{n;'DBMj 

—f - \ • 

i<Np<Mp,fxev T^^Np) 


(14) 


Note the fact that q{ii]VBmp) < ^sQih^^'B^BNp) for all 1 < Np < Mp follows due to the definition of 
Saturation Constant. Note that < 1 implies that Q{lJi]VBM ) < q{i^'i'B>Bn ) for all 1 < Np < Mp. 
In other words, p(/L 6 ; •) strictly decreases as more points are included in the database. Similarly, = 1 
implies a monotone decrease in •) as the number of points in the database increases. If > 1, 
g{pir) may increase at some point pi G VB for certain iterations of the greedy algorithm. 

Assume that the saturation constant Tq is known. At Iteration A^, let Pprofiie(Mc) denote the most 
recent available error estimate at G S: Let Q^axO^^Np) denote the maximum error estimate among 
all the previously computed error estimates at Iteration A^, that is. 


f max gtap.VBN ) for 1 < c < N'^ 
I 0 for c = 1 . 


(15) 


If TsPprofiie(Mc) < ), then g{fi^]VBN ) IS guaranteed to be less than ) due to the 

definition of the saturation constant. Thus it is not necessary to compute '^^Np) if ^profile(M c) ^ 
Np)‘ Hesthaven, et al. use this property to avoid computing error indicators at some points in 
E and save computational time in their improved greedy algorithm [24]. 

However, the saturation constant is not known a priori in general except for special cases only. For 
example, if a global ROM is constructed without truncation for a symmetric coercive elliptic problem 
and the error is measured in the intrinsic energy norm, then the saturation constant Tg is one. Hence, 
in general, the saturation assumption may not be satisfied. Nevertheless, it is still possible to use this 
concept to avoid a large number of error indicator evaluations as demonstrated below. 

In the present work, the constant Tg is one initially, then estimated and updated at each iteration of 
the greedy procedure. Based on the definition of the saturation constant (14), the following estimate for 
Tg is possible at Iteration Np: 


(2profile(Mc) 

/^c ^profile 


max 


(16) 


where ^pr^fiY(Mc) ^ preceeding available error estimate to Pprofiie(Mc) Me- The points with small 
error estimates can give a misleadingly large estimate Tg. Thus, the condition ^p^fi^(Mc) ^ ^toi is 
imposed in order to avoid those points with small error estimates in estimating Tg. 

A safety growth factor 7 is then introduced because the approximation (16) is still a lower bound 
for Tg. Algorithm 1 presents the saturation assumption-based adaptive greedy algorithm for the ROM- 
interpolation model. In Line 12 of the algorithm, Ttemp(i) is taken the maximum value among the pair 

(1, • If makes Tg > 1 throughout the greedy iterations. Similarly, as in the random greedy 

V ^profileA 

procedure, a sanity check without saturation-assumption filtering is also done after convergence of the 
greedy procedure. 

Remark. The ROM-interpolation model relies on the interpolation scheme that is used to interpolate 
ROMs in VB (see Section 5). An over-fitting issue can cause a large fluctuation for some points in V even 
though several points are added to VB. In order to prevent the over-fitting issue, it is recommended to 
use a smooth interpolation that is close to a linear interpolation. This can be accomplished, for example, 
by the multi-quadric radial basis function, which is defined in (17), with a low value of 0 as depicted in 
Figure 3. 

0(r) = \/r‘^ T 0‘^ (17) 

The performance of the saturation-based adaptive greedy Algorithm 1 is compared both with a 
classical greedy algorithm and the surrogate-based greedy algorithm in Section 6.3. 
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Figure 3: Illustration of overfitting problems and linear interpolation 


5 Interpolation of parametric ROMs 

A methodology of interpolating ROMs on matrix manifolds in [3, 1, 7] is introduced. It enables the 
computation of a reduced-order model at any point in the parameter space in real time. Section 5.1 
presents a brief overview of the consistent interpolation of the local ROMs on matrix manifolds. The 
sensitivity of the interpolation on matrix manifolds is introduced for gradient-based optimization algo¬ 
rithms in Section 5.2. Section 5.3 addresses the computation of residual-based error estimates in the 
context of ROM interpolation. 

5.1 Consistent interpolation of local ROMs on matrix manifolds 

The interpolation scheme has two major steps: 1) Identification of congruence transformations and 2) 
Interpolation on matrix manifolds. The first step is required because ROMs that are constructed for 
different parameter values (e.g., and /L 62 ) are usually not expressed in a consistent set of generalized 
coordinates [3]. In other words, the corresponding ROMs A^(/L 6 ;l) similarly and 

hr{/J^ 2 )) expressed in a consistent way so that they cannot be compared directly. Fortunately, 

congruent transformation of inconsistent ROMs can be defined by noticing that there are equivalence 
classes of ROBs for a given subspace S{/j.) = range(V(/L 4 )): 

= V(m)Q}, (18) 

where Q G is an orthogonal matrix. Considering for simplicity the case of Galerkin projection, the 

ROM defined in (9) is transformed as 

Arifi) = Q^A^(/l 6 )Q b^(M) = (19) 

Therefore, an equivalence class of ROMs can be defined as 

c(A^,b^) = {Q^A^Q, Q^b^lQ G 0{k)}, (20) 

where 0{k) denotes the set of orthogonal matrices of size k. 

Optimal transformations are then computed as follows [3]. Let V(/L 6 ;l )5 ^(^ 2)5 • • - ^ denote 

the ROBs associated with each ROM. If is a reference parameter value, then in the first step, the 
following series of orthogonal Procrustes problems are solved in order to find congruence matrices that 
transform each ROB V(/L 4 ^), c = 1 ,..., into a ROB consistent with the reference basis \'{/j.i): 

IV(Mc)Qc-V(A ii)||F 

j^ 21 ) 

subject to Q^Qc = I/c 5 c=l,...,Ap. 

The optimal solution to (21) is analytically given by 

Q; = UeZ^, (22) 
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Figure 4: Schematic representation of logarithm and exponential mappings and application to interpo¬ 
lation 


where is a singular value decomposition of (e.g., see [20]). Once the transfor¬ 

mation matrices QJ are computed, a set of rotated ROMs 

= {(Q;)^A,(mJQ*, (Q:)^b(M,)} (23) 

can be obtained via (19) for c = 1,..., Np. 

In the second step, interpolation of the consistent ROMs {A^(/L4^),b^(/L4^)}, c = l,...,A/p is per¬ 
formed on matrix manifolds leading to an approximation of {Ar{/J^'^) for a new value G V. 
Interpolatmn on matrix manifolds is necessary in order to preserve any characteristics that a ROM 
operator A^(/l 6^) might have (e.g., orthogonality, non-singularity, symmetry, or positive-definiteness). 
Interpolation on matrix manifolds can be divided into three sub-steps: 

1. logarithm mappings of the consistent ROMs to a linear tangent space 

2. interpolating mapped data in the linear tangent space 

3. exponential mappings of the interpolated quantity from the linear tangent space back to the original 
manifold. 

More specifically, let be a manifold in whose elements are characterized by properties such as 

orthogonality, non-singularity, symmetry, or positive-definiteness. Let X = Ar{/J^i) G At be a reference 
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Table 1: Exponential and logarithm mappings for the matrix manifolds of interest. 


Manifold 

xN 

Nonsingular matrices 

SPD matrices 

Logx(Y) 

Y-X 

log(YX-i) 

log(X-i/2YX-i/2) 

Expx(r) 

X + F 

exp(r)X 

Xi/2exp(r)Xi/2 


Algorithm 2 Interpolation on a manifold A4 
Input: Np matrices Yi,..., belonging to A4 
Output: Interpolated matrix Y* 

1: Set a reference point X = Yi (e.g., The interpolation process takes place in the linear space TViA^)- 
2: for c = 1,Np do 
3: Compute Fc = Logx(Yc) 

4: end for 

5: Interpolate each entry of the matrices Fc, c = 1,... ,Np, independently to obtain F* 

6: Compute Y(/L6*) = Y* = ExpY,(F*) 


point and Yc = e A4 for c = 1,... ,Np an element in the neighborhood of X. Let also F be 

an element of the tangent space 7x>M at X. The logarithm mapping Logx defines a mapping from 
an element in a neighborhood of X (i.e., A/’(X)) to an element in the tangent space TxAf. On the 
other hand, the exponential mapping Expx defines a mapping from the tangent space TxAf to A4. 
The neighborhood A/’(X) is identified by the property that for any element Y G A/’(X), the equation 
Expx(F) = Y has a unique solution F satisfying Logx(Y) = F. 

The tangent space TxAf is a vector space, so it is easier to apply any interpolation scheme in TxA^ 
than in A4. Let Fc = Logx(Yc) G TxAf for c = l,...,Np and be an interpolation 

operator that takes a set of Np distinct elements Fc, c = 1,..., A/^ as inputs and returns an interpolant 
F* G Then the exponential mapping Expx(F*) returns the element Y* G AT, which completes 

the procedure of the interpolation on matrix manifolds. The procedure of the interpolation in matrix 
manifolds can be compactly expressed in the following equation: 


Y* = = Expx \X (m*; {Logx(Y,)},^^) 


(24) 


Because the exponential mapping brings an element in back to a point in AT, the interpolation 

scheme on matrix manifolds described above produces an interpolant (e.g., Ac(/l^*)) that has the same 
properties as the interpolated points Ar{/J^c)N = (e.g., orthogonality, non-singularity, and 

positive-definiteness). Figure 4 graphically depicts the exponential and logarithm mappings in a manifold 
and interpolation in the linear tangent space. Table 1 shows the expressions of exponential and logarithm 
mappings for various manifolds of matrices such as real matrices, non-singular matrices, symmetric 
positive-definite matrices. Finally, Algorithm 2 summarizes the procedure of the interpolation on matrix 
manifolds. For more detailed descriptions of the interpolation of ROMs and RGBs on matrix manifolds, 
see [2] and [3]. 


5.2 Online computation of sensitivities 

Using analytical gradients instead of a finite difference approximations can speed up the convergence of 
gradient-based optimization algorithms. Furthermore finite difference approximations in the context of 
PDE-constrained optimization can be very expensive due to the requirement of solving the underlying 
PDE Np + 1 times if forward or backward difference is used and 2Np times if central difference is used. 
To alleviate that cost, analytical sensitivities of the interpolated ROMs are derived in this section. 

As shown in Section 2, and need to be computed to obtain gradients at a given parameter 
In the context of optimization where the high-dimensional model is replaced by a ROM as in (9) 
and (10), the sensitivities of the reduced operators and need to be computed instead. In 
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Table 2: Sensitivities of exponential mappings for the matrix manifolds of interest. 


Manifold 

Nonsingular matrices 

SPD matrices 

^Expx(r) of 

^exp(r) X" 

vl/2 Oexp(r) vl/2 

dm dm 

dm 

-/V ' o J\. ' 

dyii 


the proposed methodology, is computed by interpolation within a database of consistent ROMs as 
compactly described in Eq. (24). Note that the operators {Logx(Yc)}f^i are independent of fi. Thus, 
the derivatives of Y*(/L6) with respect to ijl become 


dY^ 

d/j. 


dfj, 


(m*) 


dExpxX 


g(^*;(L„g,(Y,))a.), 


(25) 


The derivatives of the exponential mapping associated with several matrix manifolds are provided 

in Table 2. Note that the derivatives of the matrix exponential (e.g., are required both for the 

manifolds of nonsingular matrices and SPD matrices. In order to obtain the derivatives of the matrix 
exponential one can first define a matrix 


B = 


T 

0 


dr 1 

dfii 

r 


G M 


2kx2k 


Following [30], the exponential matrix of B becomes 


exp(B) 


exp(r) 

0 


aexp(r) 

dfii 

exp(r) 


(26) 


(27) 


Hence, the derivative of the matrix exponential can be simply extracted as the (l,2)-block of 

exp(B). Note that the computation of exp(B) is inexpensive as it operates on a reduced size matrix of 
dimension 2k. 

As for and one can apply the same interpolation technique described above on the matrix 
manifold Once and are computed, the sensitivity ^ can be computed by the chain rule: 


A(Wr(Ai),Ai) = + 

djJji 


dq 

dWr 





dAr 

d/ii 


{lJ.)Wr{n) . (28) 


5.3 Computation of residual-based error indicators 

When residual-based error indicators are used in the greedy algorithm, the residual is computed as 

r(/L6, w^) = h{fjL) - A()u)V(/L6)w^, (29) 

where is the solution to A^(/L6 )w^ = b^(/L6) and Y{/j.) is the ROB. In the present context, if fJ^ ^ VB^ 
only the reduced order models (i.e., A^ and b^) are interpolated and V(/L6) is not available. It is 
possible to obtain the ROB V(/L 6 ) by interpolating the ROBs V(/L 6 ^), c = 1,..., on the Grassmannian 
manifold [2]. However, the interpolation on the Grassmannian manifold is much more expensive than 
the interpolation of the ROMs because it requires the thin SVD of large scale matrices. One way of 
avoiding the interpolation of the ROBs is to replace V(/L6) with V(/L6^) where G VB is the closest 
point to fjL (e.g., ijl^ = argmin||/L6 — The choice of is illustrated in Figure 5 for a case where the 

database VB has five points in For example, the interpolants at uses the ROB at fi^ to compute 
a residual 

r(/L6*, w^) = b(/L6*) - A(/L6*)V(/L6 i)w^. (30) 
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Figure 5: ROB choice for residual evaluation 

6 Design optimization of a wing under aeroelastic constraints 

The following design optimization problem of a wing under aeroelastic constraints is considered to 
demonstrate the performance of the proposed ROM database strategy: 

w 

D{^l) 

<^Vm(m) ^ <^upper 

W{fl) < VFupper 
Slower — M — Mupper? 

where /Lt is a design parameter vector belonging to a desing space V, L{/j.) and D{/j.) are the lift and drag, 
respectively, and W{/j.) is the weight of the wing. The von Mises stress (Jvm of the wing at the steady 
state is constrained not to exceed an yield stress crupper- The damping ratio vector C(m) is not allowed 
to be below a lower bound Ciower > 0 to avoid flutter. Bound constraints on are introduced to avoid 
unrealistic designs. Finally, two types of design parameters are considered in /j, = external 

shape parameters G and structural parameters G . The external shape parameters 
affects both the shape of the structure and the fluid domain through its interface with the structure. 
The vector of parameters fjL^ G contains the material properties of the structural system as well as 
internal shape parameters associated with the “dry” elements of the structure that are not in contact 
with the external flow. 

The constraints can be grouped into two sets: 

• Static aeroelastic constraints that are computed using a HDM that derives from a three-field 
formulation [28]. These constraints include the weight, lift-drag ratio, and von Mises stress. 

• Dynamic aeroelastic constraints (flutter) computed by a linearized HDM that derives from a three- 
field formulation linearized around an equilibrium state [27]. The proposed ROM database strategy 
is used to alleviate the large computational cost associated with the HDM for the flutter constraints. 

The damping ratio and its sensitivities ^ are then computed via the ROM interpolation technique 
described in Section 5.2. 

Section 6.1 presents the linearized aeroelastic equation and its model reduction. The section also 
shows the derivation of the aeroelastic damping ratios and their interpolations. The computational 
results regarding solutions of the optimization problem (31) are presented in Section 6.2. 


minimize 
subject to 


13 




6.1 Linearized CFD-based fluid-structure interaction and model reduction 

A CFD-based, nonlinear, high-fidelity aeroelastic system can be described by a three-field Arbitrary 
Lagrangian Eulerian (ALE) formulation that considers the moving mesh as a pseudo-structural system 
and can handle large deformations [17]. After semi-discretization in space, linearization of the semi¬ 
discrete equations about an equilibrium state [27], and the elimination of the fluid mesh position, the 
following system of linear ODEs is obtained: 


A(/l6)w + H(/l6)w + R(/l6)u + G(/L6)u = 0 (32) 

M(/L6)u + D(/L6)u + K(/L6)u = P(/L6)w. (33) 

where a dot denotes a derivative with respect to time t, A G the diagonal matrix of cell volumes 

in the fluid mesh, and Nf the dimension of the semi-discretized fluid subsystem. The conservative 
fluid state vector is denoted as w(t) G and u(t) G is the vector of structural displacements 
of dimension Ng. The Jacobians of the numerical fluxes with respect to w and x are denoted as 
H G and G G respectively. The matrix R G appears due to the linearization 

of the underlying HDM with respect to the fluid mesh velocity. Both G and R are coupling terms. The 
mass matrix M G is associated with the Einite Element (EE) discretization of the structural 

subsystem and D G ^ ^NsxNs respectively the EE damping and stiffness matrices. 

The Jacobian of the aerodynamic forces P G acts on the surface of the structure with respect 

to w. 

Eirst, the dimensionality of the structural subsystem is reduced using modal truncation. A ROB 
X(/L6) G constructed using the first kg modes of the structural subsystem and Eq. (33) 

reduced by Galerkin projection onto X(/L6). This results in an approximation of u(t) as 

u{t) ^ (34) 

where the ROB X depends on the parameters G ^ ]^ks ^g ^ vector of kg 

generalized coordinates associated with the modal approximation. The projection of (33) results in a 
set of kg equations 

ilr + D^(/L6)u^ + = X(/L6)^P(/L6)w, (35) 

where D^(/L6) = X(/L6)^D(/L6)X(/L6) G M^^x/cs ^ ^ksxks ^g diagonal matrix of eigenvalues 

associated with the kg eigenmodes. 

The dimensionality of the fluid subsystem is then reduced by POD in the frequency domain [25]. A 
ROB G is built and the state vector w approximated as 


w(t) ^ (36) 

where is a vector of kf generalized coordinates associated with V(/L6). The procedure of constructing 
V(/L6) is outlined in Algorithm 3. Note that the orthogonality condition V(/L6)^A(/L6)V(/L6) = Ikf is 
satisfied [4]. Galerkin projection of (32) using V(/L6) results in the set of reduced coupled linear ODEs: 

+ H^(/L6)w^ + R^(/L6)u^ 4-G^(/L6)u^ = 0 (37) 

ilr -\-'Dr{f^)Ur R = 0, (38) 

where H^(/l^) = V(/L6)^H(/L6g)V(/L6) and the coupling matrices are Rr(M) = V(/L6)^R(/L6g)X(/L6) and 
P^(/L6) = X(/L6)^PV(/L6). In compact form, the system (37)-(38) can be written as 


q = N(/[x)q, 


(39) 


where q(t) = 


Wr(t) 

Ur(^) 

Ur(t) 


G and 


N(ax) = 


-h,(m) 

Pr(M) 

0 


R^. {fi) 
-L)r(M) 


— Gr{fl) 

0 


G ]^(^/+2^s)x(/c/+2/cs)^ 


(40) 
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Algorithm 3 Construction of a fluid ROB for the aeroelastic system 

Input: Parameter fji = {/jL^^/jL^) G frequency sampling eigenmodes X(/L6) G R^^x/c, 

ROB dimension kf 

Output: ROB V(/L6) 

1: Compute A = A(/l6), H = R = R(/l^), and G = G{/jl) 

2: for i = 1, • • • , /cg do 
3: for / = 1, • • • , do 

4: Solve the linear system A + = — (j^^R + G)xi, where is the i-th vector in X 

5: end for 

6: end for 

7: Construct the complex-valued snapshot matrix W = [wip, • • • , 

8: Compute the singular value decomposition: A^ [Re(W), Ini(W)] = UXZ^ 

9: Retain the first kf left singular vectors V(/L6) = 


For simplicity, and without loss of generality, the case of an undamped structure is considered in the 
remainder of this paper, that is D^(/L6) = 0. For a given parameter vector fj, and altitude h, structural 
eigenvalues can be extracted from the the eigen-decomposition of N(/L 6 ) in (40): 

N(M)qj(M) = j = l,...,2ks + kf. (41) 

An algorithm for efficiently extracting the 2ks structural eigenvalues of N(/L 6 ) is introduced in [ 6 ]. It 
proceeds by extracting structural eigenvalues from a smaller size nonlinear eigenvalue problem that is 
equivalent to the eigen-decomposition of N(/L 6 ). The nonlinear eigenvalue problem is defined as 

Ns(Aj(M);M)qis(M) = 0, (42) 


where Xj is an eigenvalue and qjs the associated structural part of the eigenvector q^. The matrix 
Ng G r 2 ^sx 2 /cs defined, in turn, as 

N,(A;^) = NM - AIsfe, + (AI^^ - N//(m))“' N/,(/x), (43) 

where the blocks € ]R 2 fe.x 2 fe,^ ^ ]^ 2 k,xkf^ N// e and N/^ e defined as 


(m) 

n./(m) 


0 ■ 

. Ife. 0 J ’ 

Pr(At) 

0 J ’ 


(44) 

(45) 


(46) 

N/,(m) = [ -Rr{^l) -Gr{^l) ] . (47) 


The nonlinear eigenvalue problem (42) can be solved either by a fixed-point iterative method or a 
continuation method as described in [ 6 ]. Once a structural eigenvalue Xj is obtained, the corresponding 
damping ratio (js is defined as 



where Xf and Aj are real and imaginary part of Aj, respectively. 

In the context of solving the optimization problem (31), the operator Ns(A; fi) needs to be constructed 
and evaluated for various values of (jl. However, 'Ns{Xsi ; A^) is not available for /j, ^ V. Thus the consistent 
ROM interpolation approach described in Section 5 is applied to construct Ns(A;/L 6 ) for /j. ^ V. Let 
Nj = NgfifJ^) {Xlkf — N//(/L 6 )) ^NfsifJ^) and note that — AI 2 /C 3 + Nj by (43). Note that 

the block N 55 is symmetric positive definite, so it is interpolated on the manifold of SPD matrices. The 
matrix Nj is a singular matrix. Thus Nj is interpolated on the manifold of square matrices r 2 ^sx 2 /cs^ 
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In gradient-based optimization algorithms, the sensitivity of the damping ratio with respect to opti¬ 
mization parameter ^ needs to be provided. The detailed derivation of using consistent interpo¬ 
lation on matrix manifolds is described in Appendix 8. 


Table 3: ARW-2 specifications 


Parameter 

Value 

Geometry (inch) 

Wingspan 

104.9 

Root 

40.2 

Tip 

12.5 

Material properties 

Skin (except flaps) 

Various composite materials 

Stiffeners 

Aluminum 

E (psi) 

1.03 X 10^ 

p (Ibf-s^/in^) 

2.6 X 10-^ 

V 

0.32 



Figure 7: Variations of the shape variables 
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Table 4: Performance comparison of the various greedy algorithms considered. 


Algorithm 

Classical 

Random 

Fixed Ts 

Adaptive Tq 

Surrogate 

Ns 

125 

125 

125 

125 

125 

Nu 

- 

20 

20 

20 

10 

Number of Sampled HDMs 

16 

15 

15 

13 

15 

Max. Rel. OIBEE (%) 

4.9 

3.9 

4.6 

4.7 

4.6 

Total Elapsed Time (hour) 

60.9 

12.1 

8.8 

8.4 

7.3 

Speed-up wrt Classical Greedy 

1 

5.1 

6.9 

7.3 

8.3 



Figure 8 : Three groups of stiffeners of ARW -2 

6.2 Design optimization of the ARW2 

The optimization problem (31) is solved for the Aeroelastic Research Wing (ARW- 2 ) [33]. Physical 
dimension and material properties of the wing are reported in Table 3. A detailed FE model of the 
structure (Figure 6 ) is considered that includes, among others, spars, ribs, hinges, and control surfaces 
of the wing, and that contains a total of 2 , 556 degrees of freedom. A three dimensional unstructured 
fluid mesh (Figure 6 ) around the wet surface with 63,484 grid points is generated. An operating flight 
configuration is set for the altitude h = 4, 000 ft with atmospheric density poo = 1.0193 x 10“^ lb • s^/in^ 
and atmospheric pressure Poo = 12.7 Ib/in^, Mach number Moo = O-S, a zero angle of attack. 

The upper bound for the von Mises stress constraint (crupper) is 2.5 x 10^ psi and the upper bound for 
the weight (Wupper) 400 lbs. The lower bound for the damping ratio (Ciower) is chosen between 4.1 x 10 “^ 
and 4.5 x 10“^. Three external wing shape parameters G and three structural parameters G 
are considered. The structural shape parameters include back sweep angle (p^), twist angle (p^), and 
dihedral angle (p^), which are depicted in Figure 7. The structural parameters are thickness 

increments for three disjoint groups of stiffeners. The groups of stiffeners are depicted in Figure 8 . 

Two distinct optimization problems are considered: 1) only the external shape parameters are used 
as optimization variables (i.e., 2 ) both external shape and structural parameters are considered 

(i.e., jj, = Both cases are solved with the MATLAB fmincon implementation of the active 

set method. Both ij,^ and are also normalized so that their upper and lower bounds are 0.1 and 
—0.1, respectively. For each problem, the proposed ROM database approach for handling the flutter 
constraints is compared to an optimization based on the HDM. 

6.3 Offline database construction 

Table 4 compares the performance of several greedy algorithms in the offline phase of the ROM- 
constrained optimization when /j, = (see Section 3). The candidate set is constructed by full factorial 
design with five points each axis, leading to = 125. The subset size A^n is 20 and the size of the set 
for the sanity check 50. The marginal factor 7 is 1 and ctoi = 0.05. 
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• Classical is for the classical greedy algorithm, in which error estimates for every candidate point 
is evaluated. 

• Random denotes a greedy algorithm that chooses a random subset of the candidate set as described 
in Section 4.2.2. 

• Fixed Ts denotes a saturation-assumption based adaptive greedy algorithm where the fixed satu¬ 
ration constant = 2 is used. 

• Adaptive Tq denotes a saturation-assumption based adaptive greedy algorithm where the saturation 
constant Tq is modified at each greedy iteration according to Algorithm 1. 

• Surrogate denotes an adaptive greedy algorithm where a surrogate surface of error estimate is used 
to pick a next point and add to the database (e.g., see [31]). 

All the greedy algorithms use the One-Iteration-Based Error Estimate (OIBEE) as an error indicator 
proposed in [12]. There is a particular characteristic of evaluations of OIBEE for ARW-2 damping ratio 
when shape parameter is present. It is less expensive to evaluate OIBEE if a set of shape parameters 
is revisited because all the computations required to update the shape of ARW-2 have already been 
processed. Classical, Random, Fixed, and Adaptive greedy algorithms have higher probabilities of revis¬ 
iting a point in parameter space than Surrogate because they work on the whole or a random subset 
of a fixed set of candidate points at each greedy iteration but Surrogate do not. This is why Surrogate 
is more expensive than other three adaptive greedy algorithms in Table 4. All the greedy algorithms 
except Classical have a randomness in choosing a subset of candidate set. The results shown in Table 4 
are from one representative instance of greedy simulations. By taking a random subset of the candidate 
set, the speed-up of at least 5 is achieved from Classical. Eurther speed-up is achieved by using the 
saturation assumption filtering. A similar performance is achieved by Fixed and Adaptive Tq (e.g., a 
speed-up of 6.9 for Fixed Tq and 7.3 for Adaptive Ts). 

6.4 Speed-up for one online evaluation of the flutter constraint due to ROM 
database 

Table 5 compares the performance of the HDM and the ROM database models for one online evaluation 
of the flutter constraint {(^ and ^). Both wall-clock and CPU times are reported where CPU times are 
calculated as products of wall-clock time and the number of CPUs. Eor pL = pL^, wall time speed-up 
of 16.7 and CPU time speed-up of 535.1 are achieved. Eor pi = {pLg,pL^), wall time speed-up of 28.9 
and CPU time speed-up of 926.1 are achieved. Both cases emphasize the large computational speed-up 
achieved by the proposed ROM database strategy. 

6.5 Online predictions and optimization for N, = 3 

Table 6 presents the optimization results for pi = pi^ and Clow = 4.2 x 10“^. Both the HDM and the 
ROM database approach start from the same initial shape pi = (0.1,0.1,0.1) and converge to similar 
optimized designs. The initial and optimized shapes of ARW-2 are depicted in Eigure 9 for the HDM- 
based optimization. Although the optimized shape of ARW-2 for the ROM database approach-based 
optimization is not depicted in Eigure 9, it is almost the same as the one for the HDM-based optimization. 
The flutter constraint is violated at the initial shape, but is satisfied at the optimal shape in both cases. 
The lift-drag ratio increases by 15.1%, the weight and the maximum von Mises stresses are raised by 
4.4% and 18.2%. Although a similar number of optimization iterations is taken both for the HDM and 
the ROM database approach (5 and 6, respectively), the optimization path is different. Eigures 10 and 
II show the ARW-2 shape changes corresponding to both optimization paths. Eigure 12 shows the 
optimization iteration history of shape variable values and Eigure 13 shows the optimization iteration 
history of various quantities such as lift-drag ratio, minimum damping ratio, maximum von Mises stress, 
and weight. 

Table 6 also reports the CPU times and corresponding speed-ups. The Offline Phase reports the 
CPU hours required to construct a database of ROMs by a greedy algorithm using the adaptive Tq 
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Table 5: Computational time (hour) for one online evaluation of the flutter constraint and its sensitivities 




= Ms (-^/x = 

3) 

M = 


= 6) 


# of CPUs 

Wall Time 

CPU time 

# of CPUs 

Wall Time 

CPU time 

HDM 

32 

0.342 

10.934 

32 

0.492 

15.744 

ROM database 

1 

0.019 

0.019 

1 

0.017 

0.017 

Speed-up 


18.0 

575.5 


28.9 

926.1 


approach and the Online Static Constraint Evaluation shows the CPU hours for computing all the static 
constraints and objective functions. The Online Flutter Constraint Evaluation shows the CPU hours of 
computing flutter constraints (i.e., damping ratio and its sensitivities). A speed-up of 559.3 is gained 
when comparing the Online Flutter Constraint Evaluation. A speed-up of 17.6 is achieved when the total 
online computation is compared. This smaller speed-up is due to the fact that the static constraint HDM 
is not reduced. In the present paper, the overall total computational time for the ROM-interpolation 
approach is larger than the one for the HDM due to the expensive cost of constructing a database ROM 
in the offline phase. 

Table 6: CPU timings for an ARW-2 shape optimization problem with (^low = 4.2 x 10“^ and 


Design 

Initial 

Optimized 

Optimized 

Approach 


HDM 

ROM Database 


(0.1,0.1,0.1) 

(-0.1,-0.0814,0.1) 

(-0.1,-0.0807,0.1) 

LjD 

10.6 

12.2 

12.2 

Weight (lbs) 

342.6 

357.7 

357.7 

Min. ^ 

2.24 X 10-^ 

4.22 X 10-^ 

4.25 X 10-^ 

Max. (JvM (psi) 

18,731.9 

22,139.5 

22,139.0 

Optimization Iterations 


5 

6 

Function Evaluations 


24 

22 

CPU time (hour) 




Offline Phase 


0 

268.8 

Online Static Constraint Evaluation 


14.2 

14.1 

Online Flutter Constraint Evaluation 


240.5 

0.43 

Online Phase Total 


254.7 

14.5 

Total (Offline+Online) 


254.7 

283.3 

Speed-up 




Flutter Constraint Speed-up 


1 

559.3 

Online Phase Speed-up 


1 

17.6 

Total Speed-up 


1 

0.7 
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Sweep Angle (/ii) 


Figure 9: Initial and optimized configuration of ARW-2 


Figure 10: ARW-2 shape history of HDM optimization (Iterations 1 to 5 from the left to right) 



Figure 11: ARW-2 shape history of ROM optimization (Iterations 1, 2, 3, 4, 6 from the left to right) 



Figure 12: History of shape variables for the ARW-2 optimization 
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Figure 13: History of various quantities for the ARW-2 optimization 

However, the database of ROMs can be reused for multiple optimization problems such as multi¬ 
start, multi-objective optimization, and robust optimization. In these contexts, the ROM-interpolation 
approach has the potential to lead to much more speed-up. Table 7 presents the CPU time for multiple 
flutter optimization problems of the form (31) with nine different values of lower bounds for the damping 
ratio (Clow = {4.1,4.15,4.2,4.25,4.3,4.35,4.4,4.45,4.5} x 10“^). For each lower bound of damping ratio, 
a multi-start strategy with ten different initial points is applied to seek a global optimizer. Therefore, 
ninety independent optimization problems are solved using both the HDM and the ROM database 
approach. A speed-up of 1459.3 is achieved for the flutter constraint computations and the speed-up of 
43.6 for the online phase. A speed-up of 17.1 is achieved even for the total CPU time including the offline 
phase of constructing a database ROM and online phase. Figure 14 shows the optimal lift-drag ratio 
found from the multiple flutter optimization solutions and the corresponding weights and the maximum 
von Mises stresses. The optimal lift-drag ratio decreases as the lower bound for the damping ratio 
increases. On the other hand, the corresponding maximum von Mises stress increases. Figure 15 (left) 
shows the corresponding minimum damping ratios. The ROM database is constructed by the saturation 
assumption-based greedy algorithm with the maximum OIBEE convergence threshold of 5%. Eigure 
15 (right) shows 5% error bar and actual relative errors of the ROM database model. The minimum 
damping ratio computed by the ROM database model falls much below the convergence threshold 5%. 
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Lift/Drag ratio 


Table 7: Total CPU time for multiple flutter optimization problems 



RDM 

ROM Database 

Speed-up (hours) 

Offline 

0 

268.8 

Static Constraints 

2,065.1 

755.7 

Dynamic Constraints 

31,813.5 

21.8 

Miscellaneous Online 

14.2 

0.06 

Online Total 

33,892.8 

777.6 

Total 

33,892.8 

1,046.4 

Speed-up 

Dynamic Speed-up 

1 

1,459.3 

Online Speed-up 

1 

43.6 

Total Speed-up 

1 

32.4 



Figure 14: Optimal design dependence on the flutter constraint 




Figure 15: Minimum damping ratio of optimal design dependence on the flutter constraint 
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6.6 Online predictions and optimization for N, = 6 

Table 8 presents the optimization results for /j, = and Ciow = 4.2 x 10“^. Both the HDM and 

the ROM database models start from the same initial point and lead to the same optimal shape . The 
two models, however, lead to slightly different optimal structure material parameters (jl^. The flutter 
constraint is violated at the initial design, but is satisfied at the optimized design. The lift-drag ratio 
increases by 8% and the weight is raised by 4.7% for the ROM database model and remains almost 
constant for the HDM. The maximum von Mises stress decreases by 6.6% for the ROM database model 
and 6.4% for the HDM. Although the optimal shape parameter = (—0.1, —0.1, 0.1) is similar to the 
one in Table 6, the maximum von Mises stress for /L6 = /j,^) is smaller than the maximum von Mises 

stress corresponding to the optimal solution for fx = This is accomplished by increasing the thickness 
of the stiffeners. A similar number of optimization iterations and function evaluations are taken for the 
HDM and the ROM database model. 

Table 8 also reports the CPU time and corresponding speed-ups. A speed-up of 1154.7 is gained if 
only the dynamic constraint computation is compared. A speed-up of 14.2 is achieved if the total online 
computation is compared. The total computational time for the ROM-interpolation approach is larger 
than the one for the HDM because of the expensive cost of constructing a database ROM in the offline 
phase. However, the ROM database model can be reused in multiple optimization solves such as robust 
optimization and multi-objectives optimization problems. Figure 16 shows the total CPU time speed-up 
dependence on the number of optimization solves under the assumption that one optimization solve for 
the ROM database model (or the HDM) takes the similar time to the one in Table 8. For example, 52 
optimization solves give a total time speed-up of one. The total CPU time speed-up converges as the 
number of optimization solves increases to 14.2, which is the online time speed-up. 



Number of optimzation solves 


Figure 16: Total CPU time speed-up dependence on the number of optimization solves 


7 Conclusions 

A novel methodology for a solution of PDF-constrained optimization problems is introduced. A database 
of local ROMs is constructed and consistent ROM interpolation is performed to accelerate the optimiza¬ 
tion phase. The methodology is applied to the design optimization of a realistic aeroelastic wing. A large 
online CPU time speed-up is achieved. The online flexibility of the ROM database approach results in 
the applicability to multiple-objectives optimization, robust optimization, and multi-start strategies for 
a global optimization. In the context of multiple solutions of optimization, in turn, a total CPU time 
speed-up can also be gained. The accuracy of the ROM database model can be tuned in the offline phase 
where a database of local ROMs is constructed by a greedy procedure. A novel greedy algorithm based 
on the saturation assumption, in which a value of the saturation constant is adaptively updated in each 
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Table 8: CPU timings for the ARW-2 shape optimization problem with ^low = 4.2 x 10 ^ and = 


Design 

Initial 

Optimized 

Optimized 

Approach 


HDM 

ROM Database 


(0,0,0) 

(-0.1,-0.1,0.1) 

(-0.I,-0.I,0.I) 

Mm 

(0,0,0) 

(0.1,0.1,-0.015) 

(0.1,0.1,0.013) 

L/D 

11.3 

12.2 

12.2 

Weight (lbs) 

349.9 

349.8 

366.3 

Min. 

3.7 X 10-^ 

4.2 X 10"^ 

4.2 X 10-^ 

Max. (JvM (psi) 

20,297.1 

18,991.9 

18,953.3 

Optimization Iterations 


6 

5 

Function Evaluations 


II 

9 

CPU time (hour) 




Offline Phase 


0 

8,800 

Online Static Constraint Evaluation 


II.O 

12.8 

Online Elutter Constraint Evaluation 


173.2 

0.15 

Online Phase Total 


184.2 

13.0 

Total (OfflineTOnline) 


184.2 

8,813 

Speed-up 




Elutter Constraint Speed-up 


I 

II54.7 

Online Phase Speed-up 


I 

14.2 

Total Speed-up 


I 

0.02 


iteration, is also introduced and shown competitive performance compared to other greedy algorithms. 


8 Appendix 

Differentiating (42) with respect to /ij, z = 1, • • • , Nf^ leads to 

ONs , dNs dXj dqjs 

“q—Q js + ox o—~ 

dfii dXj dni dni 

Multiplying by the left eigenvector from the left to (49) and noting that = 0 leads to 

5N,, 


(49) 


dni 


Pfa, 




-qja 




aN, 

dX. 


(50) 


-qjs 


The sensitivity of the damping ratio Qjg with respect to fii is then obtained by chain rule: 

^ ^ ^ f jXfr _ J_\ dX^j XfX] 
dUi djii |Aj|3 \Xj\j diJi\Xj\^' 

Taking the derivative of (43) with respect to fii leads to 


(51) 


( 52 ) 
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dNf 

The sensitivity ^ ^ of the operator Ny interpolated on the matrix manifold l^g obtained 

by (25) and Tables 1 and 2: 


dNf _ dl 
d^ii 




(53) 


The sensitivity 


ON,, 

diii 

SN 


of the operator interpolated on the manifold of SPD matrices is 
aexp (X (m*; {log(N-y 2 N,,^N-y 2 }f^i)) 


dfii 


= Nyy 


djii 


-Njy, 


(54) 


where the derivative of the exponential matrix can be obtained by (27). Taking the derivative of (43) 
with respect to \j leads to 

) I. 


dXj 


(55) 


Since Ny(A,-,^*) = Ny(A,-,^i) +x(/x*;{N/(A,-,^J -N/(A,-,Mi)}^^), 

^(A,-,Ai*) = ^(A,-,Mi) + ^(M*;{N/(A,-,Mc)-Ny(A,,A.i)}^^) 




(56) 


dXi 


ic=r. 


The second equality in (56) is true if the interpolation operator X [jjX] {Tc}^!^ is linear in Tc, c = 

1,..., A/p. Indeed, this property is true for many interpolation operators, including polynomial and RBF 
(Radial Basis Function)-based interpolations. Note also that Nj(Aj, = N5j(/L4^)(AjI—N jj(/L4^))“^N fs{/^c) 

for c = 1, • • • , Np. Hence 


dNf 




(57) 


The derivative is then obtained by plugging (57) in (55). 
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